Fracture precursors in disordered systems 
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A two-dimensional lattice model with bond disorder is used to investigate the fracture behaviour 
under stress-controlled conditions. Although the cumulative energy of precursors does not diverge 
at the critical point, its derivative with respect to the control parameter (reduced stress) exhibits 
a singular behaviour. Our results are nevertheless compatible with previous experimental findings, 
if one restricts the comparison to the (limited) range accessible in the experiment. A power-law 
avalanche distribution is also found with an exponent close to the experimental values. 
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Fractures are very complex phenomena which involve 
a wide range of spatial and sometimes temporal scales. 
Accordingly, the development of a general theory is quite 
an ambitious goal, since it is not even clear whether a 
continuous coarse-grained description makes sense; ad- 
ditionally, for the very same reason, realistic simulations 
are almost unfeasible. However, such difficulties have not 
prevented making progress on several aspects of fracture 
dynamics such as propagation velocity, roughness, or the 
failure time under a constant stress 0,12,0]. In this paper 
we are interested in studying the development of the so- 
called precursors, microcracks preceding the macroscopic 
fracture in a brittle disordered environment. Some recent 
experiments [H 0- El suggest that we are in the presence 
of a critical phenomenon, although the accuracy is not 
yet high enough not only to discuss its universality prop- 
erties, but also to assess the order of the transition. 

Fractures are typically studied either by increasing 
the strain or the stress and recording the acoustic emis- 
sions generated by the microfractures preceding the final 
breakup. In both cases, there is experimental evidence 
that the probability density N(e) of microfractures with 
energy between e and e + de, follows a power law 

N( E )~e-P, (1) 

although it has been found that different materials are 
characterized by different values of the exponent (3: 1.3 
in synthetic plaster 0,1-51 in wood 4], 1.9 in fiberglass 
and 1.25 in paper |J. 
On the other hand, the cumulative energy emitted 
while approaching the fracture exhibits a qualitatively 
different behaviour depending whether the strain or the 
stress is controlled. In the former case, after a bond 
breaks, lattice rearrangements lead to a stress reduction 
which, in turn, increases the overall stability. As a conse- 
quence, a critical behaviour can be realistically observed 
only in the latter case. For instance, in Refs. 0,0 it was 
found that 

E{Pr) ~ {^pP) X ■= Pr" , (2) 




FIG. 1: Sketch of the triangular lattice. The solid lines refer 
to the nearest neighbour interactions involving the site de- 
noted by a full triangle. Left and right boundaries are pulled 
apart by a force ±<ru. Periodic boundary conditions are as- 
sumed along upper and lower borders (e.g., the point A can 
be identified with B). 



where E is the cumulative energy released up to pressure 
P, P c is the critical pressure corresponding to the macro- 
scopic failure, and P r is the so-called reduced parameter. 
In the past, simplified models have been introduced in the 
hope to capture the essential ingredients of the process. 
Among them, random fuse networks have been quite pop- 
ular, since only one scalar variable (the electric current) 
is needed to describe the onset of a macroscopic failure @. 
After the concept of self-organized criticality (SOC) was 
introduced by P. Bak the possible interpretation of 
fractures as instances of such critical phenomena became 
an appealing perspective to many researchers. As a re- 
sult, probabilitstic models inspired by the simplest SOC 
ideas have been introduced and numerically investigated 
[Io| . Altogether, it has emerged a scenario of the fracture 
as a critical phenomenon, but the relative "distance" be- 
tween models and physical reality leaves doubts about 
the validity of such conclusions. 
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More realistic models where each particle feels the force 
held induced by its nearest neighbours were already intro- 
duced in the '80s[H[H[l![I|, 

when preliminary stud- 
ies of small 2-dimensional lattices with quenched disorder 
have been carried out. However, only in the late '90s it 
has become possible to simulate sufficiently large systems 
to attempt a scaling analysis. In Ref. 0] a square lattice 
with two- and three-body interactions was studied, find- 
ing (3 = 2.5 in agreement with the mean-field behaviour 
of the fiber bundle model 0, 0] , but larger than the 
experimental values. On the other hand, large statistical 
uncertainties prevent drawing any conclusion about the 
scaling behaviour of E(P r ). 

In this Letter, studying a slightly simpler model [13, 
we are able to investigate the behaviour of E{P r ) obtain- 
ing results that are compatible with the Lyon experiment 
4J. However, having here investigated a wider range of 
parameter values, we are led to exclude a divergence of 
E(P r ) near the critical point. 

Our model consists of point-like particles sitting on 
a 2d triangular lattice with nearest-neigbour interactions 
mediated by central forces (see solid lines in Fig^l. More 
precisely, the force acting on the ith particle due to 
the interaction with the jth one is 

fy^dr.-r.D^L , (3) 

where denotes the ith particle position, | • | repre- 
sents the modulus operation. Moreover, Fij(x — a) — 
0, with a being the lattice spacing at equilibrium and 
kij = —dFij/dx\ a > so as to ensure stability. More- 
over, as soon as the mutual distance |r^ — r ^ | becomes 
larger than some threshold a*j, the interaction strength 
is irreversibly set equal to 0. As the a*-'s are gener- 
ally chosen close to a (see also Ref. 0), only small de- 
viations from equilibrium can be expected and, accord- 
ingly, the force can be linearized. Here, we have however 
preferred the equivalent choice (at the first order level) 
F(x) — — (x 2 — a 2 )/(2a), since it avoids the computa- 
tion of a square root for each bond, thus ensuring faster 
simulations. 

Furthermore, an external force au (— cru) is applied to 
each particle of the right (left) boundary along a direc- 
tion u orthogonal to the edges, while periodic boundary 
conditions are assumed along the upper and lower bound- 
aries. Since we are interested in investigating the fracture 
process in a nearly static regime, a damping term —71", 
has been added to the forces acting on the ith particle. 
The last key element of the model is disorder, which can 
be introduced by assuming that either the bond strengths 
or the thresholds a*j are distributed according to some 
law. Here, we assume a*j = a* and a dichotomic distribu- 
tion of bond strengths: a fraction c of them is set equal to 
from the very beginning, while the rest are set, without 
loss of generality, equal to 1 . As pointed out in Ref. [13J , 
this choice is equivalent to a dichotomic distribution of 
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FIG. 2: The stress a as a function of the strain d, for one 
disorder realization in a lattice with 40 x 42 particles. Part 
of the curve is zoomed in the inset so as to identify single 
avalanches (the numbers over the horizontal steps represent 
the number of broken bonds in each avalanche). 



thresholds. 

It is convenient to express all variables in adimensional 
units, since this helps scaling out some parameters. In 
particular, the lattice spacing a and the particle masses 
can be both set equal to 1. This implies that the modulus 
a of the applied force can be identified with the stress. 
A parameter that cannot be scaled out is the threshold, 
here fixed equal to a* = 1.04. This choice, besides being 
compatible with the linearization of the force field, allows 
us reaching sizeable large-scale deviations from a purely 
crystalline structure (we have been able to study lattice 
sizes up to L = 80). 

The equations of motion have been integrated by us- 
ing a Runge-Kutta algorithm with 7 = 0.7 and a time 
step equal to 0.3. We have verified that this choice 
guarantees the fastest convergence (in CPU-time units) 
to the asymptotic state 19] . Finally, we have chosen 
to fix the fraction of initially missing bonds equal to 
c = 0.3, a value close to, but definitely below, the rigidity- 
percolation threshold (c r = 0.3398). 

Numerical experiments consist in monitoring several 
observables while the stress is increased until a macro- 
scopic failure occurs at a — a c . In order to ensure that 
the sampled configurations remain close to equilibrium 
during the whole stretching process, (i) a is slowly in- 
creased, (ii) additional relaxation loops are allowed when 
new bonds break. One o f the relevant observables is the 
strain d = |(rj — r?) • u|, where the over line denotes the 
average over all particles of the left and right edges and 
r° represents the equilibrium position. As a typical ex- 
ample of the observed phenomenology, in Fig. we plot 
the stress-strain curve for one realization of the disor- 
der. After an initially monotonous growth, a series of 
steps follows, each corresponding to the "simultaneous" 
breaking of s bonds. Since the same amount of energy is 
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FIG. 3: Average cumulative fraction of broken bonds p versus 
the reduced stress </> for L — 20 (triangles), 40 (plusses), and 
80 (circles). The dashed curve corresponds to the scaling 
behaviour observed in Ref. Q. In the inset, 1 — p is plotted 
versus (f> with the same meaning of the symbols. 




FIG. 4: Probability density P(s) of avalanches preceding the 
macroscopic failure, scaled to its maximum value. Dashed, 
dashed-dotted, and solid lines refer to L = 20, 40, and 80, re- 
spectively. The straight line, corresponding to a decay s -1 ' 75 
is the fit - shifted for clarity - of data for L = 80 



released in each bond breaking, it is natural to identify s 
with e and interpret the phenomenon as the occurrence 
of an avalanche. 

Another relevant observable is the number N(a) of 
bonds broken at stresses smaller or equal to <r, a quantity 
analogous to the cumulative energy E in Eq. @. More- 
over, having the possibility to average over many differ- 
ent realizations of the disorder (actually, we have studied 
1000 realizations for each lattice size), it is also neces- 
sary to define a meaningful way of performing ensemble 
averages. We have decided to consider the fraction 



P 



N(t) 
N{a c ) 



(4) 



as the independent variable (N(a c ) being the total num- 
ber of broken bonds in each realization) and to average 
the stress values where p is attained. Moreover, analo- 
gously to Ref. , where the reduced pressure P r has been 
introduced, here we define the average reduced stress 



(5) 



where (•) denotes an average over disorder. Our simula- 
tions indicate that, in the thermodynamic limit L — > oo, 
the average becomes irrelevant. Indeed, we have verified 
that a c is a self-averaging quantity, by observing that the 
normalized variance 
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(6) 



decreases as A^ ~ L 2 / 3 with increasing L. 

The resulting behavior for L = 20, 40, and 80 is plot- 
ted in Fig. El The data for the L = 40 and 80 nicely 



overlap, indicating that finite-size effects are already neg- 
ligible for L ps 40. In comparison with the experiment 
of Ref. J4J, the reduced stress here covers a three-times 
wider range; this allows us ruling out the existence of 
a power-law behaviour of p. On the other hand, if we 
restrict the analysis to the range accessible in the exper- 
iment, we do not find relevant differences with respect to 
the experiment itself, as indicated by the dashed line in 
Fig. |3| Next, one notices that p remains strictly smaller 
than 1 even close to <j> = 0, suggesting that the frac- 
tion of bonds broken in the macroscopic failure is finite. 
However, the most refined presentation in the inset re- 
veals that, though slowly, p eventually converges to 1. 
Indeed, by fitting the small-0 region with p — po — bcfp , 
we find that po increases with L and is equal already to 
0.94 for L — 80. As for the rj value, our best estimate is 
r/ = 0.12±0.0I. Although the fit is quite good, the small- 
ness of 77 leaves doubts about the effective behaviour of 
p in the vicinity of the macroscopic failure. 

The dependence of p on sheds light on the critical 
behaviour around the onset of the macroscpic fracture. 
It is interesting to investigate also the dependence of the 
total number of broken bonds N(a c ) on L. Here, we 
have studied too a few sizes, to convincingly determine 
the scaling behaviour; however, assuming that eventually 
N(a c ,L) = L a , we find that a w 1.5, a value that is 
compatible with the simulations of the same system with 
imposed strain |l8( . 

Another characterististic of fracture processes often 
studied both in experiments and simulations is the distri- 
bution of microfractures (i.e. avalanches). In the present 
context, this amounts to computing the number N s of 
avalanches of size s; in order to increase the statistics, we 
sum over all disorder realizations. The probability distri- 
butions for L — 20, 40, and 80 are plotted in Fig.QJ where 
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FIG. 5: The probability density Q(y) of the final avalanche 
size y. The size is shifted around the mean value and scaled 
to ensure a unit variance. Crosses, triangles, squares and 
circles correspond to L — 10, 20, 40, and 80. The solid line 
corresponds to the Gumbel distribution. 

N s is scaled to its maximum value (P(s) = N s /N™ ax ). 
All curves exhibit a seemingly power-law decay followed 
by a peak at large s, which is clearly a finite-size effect, 
since it corresponds to avalanche sizes that are compa- 
rable with the lattice size. By fitting P(s) in the mean- 
ingful s range with as~ 13 , we find f3 = 1.75 ± 0.06, a 
value that lies inside the interval corresponding to the 
various experimental measures, [1.25,1.9], [2(J- Finally, 
we have investigated the size distribution of the macro- 
scopic avalanche. The results plotted in Fig. reveal 
increasing deviations from a Gaussian behaviour, when 
L is increased. The Gumbel distribution, derived to de- 
scribe extreme- value statistics (see the solid line) appears 
to provide a more convincing description of the avalanche 
distribution. 

Once the fracture process is recognized to be a non- 
equilibrium phase transition, it becomes important to 
establish the order of the phenomenon. From the ex- 
isting literature, it is unclear whether we are in front of a 
second order transition or the spinodal point of a first or- 
der transition [l5ll2l|. Such an uncertainty is mainly due 
to the fact that a convincing order parameter has not yet 
been identified. On the basis of our simulations and of 
the Lyon experiment, it seems reasonable to consider the 
cumulative fraction of broken bonds as a proper order pa- 
rameter. As a consequence, since the contribution of the 
final avalanche is increasingly negligible in the large size 
limit, we are led to conclude that the fracture is indeed 
a continuous transition. However, it seems necessary to 
investigate yet larger systems in order to convingly iden- 
tify the asymptotic scaling behaviour of p in the vicinity 
of the critical point. 

While we are quite confident about the validity of our 
conclusions in this specific model (since the same scenario 
arises when both the lattice geometry and the direction 



of the applied stress is changed), we should add that a 
different behaviour has been found for rectangular distri- 
bution of elastic constants. In such a case, although the 
numerical results are again compatible with the experi- 
mental ones (within the experimentally accessible range) , 
p(<p) seems to exhibit a discontinuity at <fi = 0. It thus not 
illogical to conjecture that the fracture may fall within 
(a few) different universality classes, depending on some 
microscopic details of the model that have still to be iden- 
tifed. 
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